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The structure of a uniform simple liquid is related to that of a reference fluid with purely repulsive 
intermolecular forces in a self-consistently determined external reference field (ERF) (pR. The ERF 
can be separated into a harshly repulsive part ^jjo generated by the repulsive core of a reference 
particle fixed at the origin and a more slowly varying part (;/>_ri arising from a mean field treatment 
of the attractive forces. We use a generalized linear response method to calculate the reference fiuid 
structure, first determining the response to the smoother part 0_ri of the ERF alone, followed by the 
response to the harshly repulsive part. Both steps can be carried out very accurately, as confirmed by 
MD simulations, and good agreement with the structure of the full LJ fluid is found. 
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I. INTRODUCTION 

In this paper we describe a new and physically moti- 
vated way to determine the structure of uniform fluids, 
based on a mean field treatment of the attractive inter- 
molecular interactions. To apply this approximation to 
a general nonuniform system, attractive interactions are 
replaced by a spatially varying single particle "molecular 
field" potential, chosen to take account of variations in the 
average attractive energy density in different parts of the 
system Since the attractive interactions usually op- 
erate over an extended range, it seems plausible that an 
averaged description as given by mean field theory could 
often provide a useful simplification. 

However, in real liquids, additional very important "ex- 
cluded volume" correlations are generated by the short- 
ranged and harshly repulsive intermolecular forces [^|J^]; 
these cannot be accurately described using the same mean 
field averaging appropriate for the longer ranged attrac- 
tive forces. Despite this additional complexity, the use of 
mean field theory allows us to consider an inherently sim- 
pler system: a nonuniform reference fluid. This consists 
of particles interacting only with repulsive intermolecular 
forces but in the presence of an external reference field 
(ERF) chosen self-consistently to take account of the lo- 
cally averaged effects of attractive interactions as well as 
any imposed external field |^-^. The uniform reference 
fluid is stable over the entire range of densities from va- 
por to liquid, and its structure in the presence of an ap- 
propriately chosen ERF approximates that of the original 
system. 

In previous work we showed how these ideas can 
be used to give an accurate description of the structure 
of a nonuniform Lennard-Jones (LJ) fluid in a number of 
different applications where the ERF is large and attrac- 
tive forces strongly influence the structure, including the 
liquid-vapor interface and the structure of fluids near hard 
walls. In these cases conventional (singlet) integral equa- 
tion methods have given poor results . 



In this paper we consider a different limit, that of the 
uniform LJ fluid. Here attractive forces produce relatively 
small structural changes at high density and integral equa- 
tion methods have had their greatest successes. Indeed, in 
the simplest picture, the attractive forces on a given par- 
ticle from oppositely situated neighbors essentially cancel 
[H in typical high density configurations, and the structure 
of the dense uniform LJ fluid is rather well approximated 
by that of a uniform reference fluid at the same density 
The theoretical challenge is to improve on this rather 
accurate starting point at high density and to describe the 
larger structural changes attractive interactions induce at 
lower densities. 

From this perspective the uniform LJ fluid provides an 
important and nontrivial test of our general approach. It 
is not clear that mean field averaging along with the ap- 
proximate methods j6|j^ we use to calculate reference fluid 
structure will be accurate enough to determine the small 
but subtle changes induced by attractive interactions in 
the highly oscillatory structure of uniform fluids at high 
density or the more substantial changes seen at lower den- 
sities. Indeed, unlike the previous applications, it is diffi- 
cult to guess even qualitative features of the ERF. 

The plan of the paper is the following. In Sec. II we 
define the nonuniform reference fluid and the formal equa- 
tion determining the ERF. In Sec. Ill we discuss the usual 
mean field approximation for the ERF and suggest a new 
generalized equation that incorporates exact results at low 
density. 

In Sec. IV we first use computer simulations to carry out 
the determination of structure essentially exactly. This 
allows us to test the accuracy of the basic mean field de- 
scription of the ERF without any further approximations. 
We generally find quite satisfactory results though small 
errors in the simplest mean field determination of the ERF 
can be seen at high density. 

In Sec. V we introduce a new theory to calculate self- 
consistently both the ERF and the associated structure in 
the reference fluid, using a generalization of a physically 
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motivated procedure first used to calculate the structure 
of a L J fluid near a hard wall . The key idea is to divide 
the ERF into rapidly and slowly varying parts. We de- 
termine the response of the reference fluid density to each 
component of the ERF in successive steps, using appropri- 
ate methods in each step that can accurately describe the 



very different density responses, as discussed in Sees. VB 



and V C . Sec. Vl discusses the results of the method and 
comparison to simulations. The theory generally gives 
quite satisfactory results. However, at the highest densi- 
ties some small errors can be seen arising from the simple 
mean field treatment of attractive interactions. At lower 
densities, our simplest approximations give results compa- 
rable to the best integral equ ation methods Final 
remarks are given in Sec. VII. Some technical details are 
presented in the Appendices. 



full system reference system repulsive system 
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FIG. 1. Model systems considered, a) is the full LJ sys- 
tem with a LJ particle fixed at the origin as indicated by the 
dashed circle. The interaction with the other particles can be 
described by an external field 4>Lj{r) = w{r). b) is the nonuni- 
form reference system with the special wall particle fixed at the 
origin with interaction cpnir). c) is the original WCA repulsive 
force system with pair potential uo(r). Here the fixed particle 
interacts with the other particles through (f)Q{r) — uo{r). 



II. NONUNIFORM REFERENCE FLUID 

We first consider the case where fiuid particles interact 
with a known external field (f)(r) and through the LJ pair 
potential w{r) = uoir) + ui{r), separated as usual |^ into 
rapidly and slowly varying parts so that all the repulsive 
intermolecular forces arise from uq and all the attractive 
forces from ui We assume that the external field 

(j>{r) = 00 (r) + 01 (r) can be separated in a similar way, 
where the subscript denotes a harshly repulsive inter- 
action and the subscript 1 a much more slowly varying 
interaction usually associated with attractive forces. We 
consider a grand ensemble with fixed chemical potential 
, which determines , the uniform fluid density when 
= 0. 

We relate the structure of the nonuniform LJ system 
to that of a simpler nonuniform reference fluid jj-]^], 
with only repulsive intermolecular pair interactions uq (rij ) 
(equal to the LJ repulsions) in a different effective refer- 
ence field (ERF) (f>ii{r). The replacement of attractive pair 
interactions by an approximate local "molecular field" is 
an essential step in mean field theory, but we can think 
of other more general prescriptions for 0/{(r). Here we 
determine 0fi(r) formally by the requirement that it has 
a functional form such that the local (singlet) density at 
every point r in the reference fluid equals that of the full 
LJ fiuid §: 

po(r;[0fl]) =p(r;[0]). (1) 

The subscript in Eq. (|^) reminds us that the reference 
system pair interactions arise only from uq and the nota- 
tion [(/jr] indicates that all distribution functions are func- 
tional of the appropriate external field. Since p(r; [0]) is a 
physically realizable distribution function, and the refer- 
ence fluid is stable over a wide range of densities, it seems 
very plausible that such a choice for the held 0^ can be 
made in principle p^ . In practice we will make approxi- 
mate choices motivated by mean field ideas. 



Using this perspective, let us consider the problem of 
determining the effects of attractive forces on the structure 
of uniform fluids. The simplest (WCA) approximation 
[|j assumes complete cancellation of attractive forces and 
approximates the radial distribution function g(r) of the 
uniform L J fiuid by the go (r) of the uniform repulsive force 
fluid at the same density. To improve on this, we make use 
of the exact relation jl^] between gir) in the uniform 
LJ fluid and the singlet density in a nonuniform fluid with 
a particle fixed at Tq, which we take to be the origin of our 
coordinate system: 



P'^giri) = p(ri|ro; [0 = 0]) = p(ri; [0lj]) 



(2) 



Here p{vi\vq] [0 = 0]) is the conditional singlet density — 
the density in zero external field at ri given that a particle 
is fixed at rg. By symmetry this depends only on the radial 
distance ri = |ri| from the fixed "wall particle" at Tq = 0. 
This in turn must equal the nonuniform singlet density 
induced by the external field 4)lj{ti) = w{ri). 

By choosing 0i?(ri) in Eq. (^ to fit the nonuniform 
LJ density p(ri; [0lj]), we obtain a nonuniform reference 
system in which the density poiri; [4>r\) is modified by the 
effects of attractive forces. In particular this can be used in 
Eq. (H) to calculate the radial distribution function giji) of 
the uniform LJ system. The original WCA approximation 
[D arises from the particular choice 0^ — uq. See Fig. (1). 



III. MEAN FIELD APPROXIMATION FOR THE 

ERF 



In previous work ^-g] we started from the balance of 
forces as described by the exact Yvon-Born-Green hierar- 
chy Q and arrived at a generalized mean field equation for 
the ERF by a series of physically motivated approxima- 
tions. We will not repeat these arguments here and instead 
focus on the simplest final approximation, the molecular 
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field equation for the ERF, which proved surprisingly ac- 
curate in a number of different apphcations. This is just 
a transcription of the usual molecular field equation for 
the Ising model |l[ to a continuum fluid with attractive 
interactions ui{r) and can be immediately written down: 

^nn) = <l>Lj{ri) + [ dr2 [po(r2; [4'^]) - P^] u,{n,) . 



(3) 

The effective field fp^^ at a particular distance ri from the 
fixed wall particle is comprised of the bare field (j)Lj{ri) 
from the fixed particle plus the integral over all positions 
T2 of the attractive interactions Ui(ri2) from other parti- 
cles weighted by the deviation of the nonuniform reference 
density po(^2; from its limiting value . Use of the 

density deviation ensures that ipR^ vanishes at large ri. 
Let (j)s denote the second term on the right in Eq. (^: 



(/)s(ri) = / dr2 [po(f2; 



(4) 



It provides an estimate of the averaged effects of attractive 
pair interactions arising from the other (mobile) particles 
in the full LJ fluid at a distance ri from a particle fixed 
at the origin. Because of the convolution with the slowly 
varying attractive potential "weighting function" ui(ri2) 
in Eq. (Q), 0s (''1) extends smoothly into the repulsive core 
region of the wall particle where po{ri; [4'r^]) vanishes. 
Outside the core it is a smooth, basically repulsive and rel- 
atively slowly varying interaction even when po{ri; [4>^^]) 
itself has pronounced oscillations. 

More complicated, but sometimes more accurate, equa- 
tions for the ERF are available H, but in practice the 
simple mean field approximation (|3|) often gives quite sat- 
isfactory results. In Appendix A we discuss a simple mod- 
ification of Eq. that gives somewhat more accurate 
results at low density. In the following we will use Eq. 
to determine the ERF unless otherwise indicated. 



IV. RESULTS FROM MD SIMULATIONS 

We now must solve Eq. (||) to determine the ERF cj)^^ 
and associated density po{r: [0^^^]). As is typical in mean 
field theory, a self-consistent solution must be found, since 
the ERF (f)ji appears explicitly on the left side and implic- 
itly on the right side through the dependence of the density 
Po(7'2; [0-r]) on If we can find the reference structure 
Po{r] [(j)fi\) produced by a given ERF (pn accurately, then it 
is straightforward to solve the mean field equation (by it- 
eration, for example) to determine the self-consistent (f)^^ 
and the associated density po{r; In Sec. ^ we will 

discuss new theoretical methods to calculate pQ^r; [(j)R\) 
for a given (f)f(_. However, since these could introduce ad- 
ditional errors, it is useful first to assess the accuracy of 



the basic mean field equation (||) without any further ap- 
proximations. 

To that end, we carried out MD simulations for the 
three model systems shown in Fig. (1): the full LJ sys- 
tem, the WCA repulsive force system, and the inhomo- 
geneous reference system with the special wall particle 
fixed at the origin [|l5|. To determine the effective po- 
tential in the latter case, we solved Eq. (^) by iteration 
[ p!6| , using (essentially exact) MD results for po{r; [0h]). 
The errors in po{r; [(pji^]) when compared to p{r; [cpLj]) 
then arise solely from the mean field approximation for 
the ERF (j)^^. 

A. Simulation details 

In the following we use reduced Lennard-Jones units 
where the unit of length is cr, the unit of energy is e and 
the unit of time is ^ma^ je. We carried out MD simu- 
lations in the (NVT)-ensemble using the velocity Verlet 
algorithm with a time step of At = 0.001. To maintain 
constant temperature, every 150 MD steps we chose new 
velocities for all particles from the corresponding Boltz- 
mann distribution. 

We simulated states along the near critical isotherm at 
T = 1.35 for densities p^ = 0.78, 0.54, 0.45 and 0.1, and a 
state near the triple point with T = 0.88 and p^ — 0.85. 
We used N = 3000 particles for T = 1.35, p-^ = 0.78 and 
for T = 0.88, p^ = 0.85, and N ^ 450 for all other states. 
To eliminate the possibility of finite size effects we made 
test runs for T = 1.35 and p^ = 0.54,0.45 and 0.1 with 
N = 3000 particles, which led to the same density distri- 
butions as for N — 450. Each state was first equilibrated 
for 5 • 10^ MD steps. Subsequently we calculated g(r) for 
the uniform systems and po{r; [<j>ji^]) for the nonuniform 
reference system for at least 3.5 ■ 10^ and up to 7.5 • 10^ 
MD steps. 



B. Simulation results for the ERF 

We first concentrate on the high density state p^ — 0.78 
and T = 1.35, which will illustrate many basic features of 
the mean field approach. Figure (2) gives simulation re- 
sults for the full LJ density p^g{r), the WCA reference 
density p^go{r), and the inhomogeneous reference density 
po{r; [4>R^]) for this state. We see that the mean field 
prediction po{r; [(t>]{^]) ~ P^9{t) is able to correct the 
main quantitative errors in the already rather accurate 
WCA approximation goir) ~ g{^)^ describing in partic- 
ular the slight shift outward of the first peak. However, 
some small errors remain, due solely to the mean field ap- 
proximation for the ERF from Eq. (^ . These are focused 
on in the inset to Fig. (2), which compares the density 
change po{r] [(t>'j(^]) — p^go{r) due to attractive forces as 
predicted by Eq. (ph to the actual change p^[g{r) — go{r)] 
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FIG. 2. Densities with a particle fixed at the origin as de- 
termined by MD simulations for the three systems in Fig. 1 
at T = 1.35 and = 0.78. The inset gives the difference in 
density between the LJ fluid and the WCA repulsive force fluid 
as determined by simulations and by self-consistent solution of 
the mean field equation (|3|), again using simulation data. 



given by the simulations. Any further improvements in 
these results will require a better approximation for the 
ERF. 

Figure (3) shows the corresponding self-consistent ERF 
(j>^^{r) from Eq. (||), compared to the bare LJ potential 
w{r) and the repulsive reference potential uo(r). At low 
density (j)^^ reduces exactly to w, and if the force can- 
cellation argument were exact, then at high density 4>^^ 
would equal uq as assumed in the WCA approximation. 
However, there is a weak negative region in (t>^^{r) for r 
between about 1.1 and 1.4. This results from the nonuni- 
form attractive energy density experienced by a particle in 
this region in the LJ system, which is slightly lower than 
average because of the fixed particle and its neighbors even 
at this high density. This produces the slight shift in the 
first peak noted above. 

Figure (4) shows the ERF for a series of states along 
the T ~ 1.35 isotherm. The attractive force cancella- 
tion from further neighbors becomes increasingly less effec- 
tive at lower densities, and attractive interactions produce 
much larger structural changes, as will be shown below. 

Again mean field theory can yield accurate results. This 
is illustrated in Fig. (5) for the low density state = 0.1, 
and T — 1.35. This figure shows simulation results for 
the full LJ density p^g{r), the WCA repulsive fluid den- 
sity p^go{r), and the inhomogeneous reference density 
Polr; [(t)R^])- Also shown is poir; [0^ ^^]), with the ERF 
calculated from Eq. ( |Al[ ) in Appendix A (using l2{p) as 
the interpolation function), which does a slightly better 
job at reproducing the second peak than does Eq. (||). 
Theoretical values for correlation functions for all the 




FIG. 3. The self-consistent ERF (j^R" {r) determined from 
the mean field equation (^) for T — 1.35 and = 0.78 com- 
pared to the full LJ potential io(r) and the repulsive force ref- 
erence potential uo{r). 



states in Fig. (4) and comparison to results for the full 
LJ fluid will be discussed in Sec. 



Vl below. 



As the density tends to zero, the ERF reduces to the 
bare potential w{r) as correctly predicted by Eq. (^). The 
increased "screening" of attractive forces as the density is 
increased was first demonstrated using diagrammatic re- 
summation techniques in the derivation of the ORPA and 
EXP integral equations ||l^ . Mean field theory provides a 
very simple and physically suggestive way of understand- 
ing these results. 



V. TWO STEP METHOD 



We now discuss theoretical methods [|6[|9| for determin- 
ing the density po(r; [0_r]) = Pfl(r) produced by a given 
ERF We use a generalization of the two step method 
first introduced in Ref. [|| . Initially we treat the L J repul- 
sive potential wg as a hard core interaction with diameter 
d, but then use the standard "blip function" expansion 
to correct for its finite softness in our final numerical 
results. Let us concentrate on the mean field equation (^) 
for the ERF. Recall that 4>Lj{r) = w(r) = uo{r) + uiir). 



A. Separation of ERF 

The two step method introduces a similar division of 
the full ERF, 



(5) 



and determines the density response to each part of the 
ERF in separate steps. As this notation suggests, 4>m is 
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FIG. 4. The self-consistent ERF from Eq. (||) for the indi- 
cated densities along the T = 1.35 isotherm. 

supposed to take account of the harshly repulsive part of 
the ERF associated mainly with the repulsive core of the 
fixed wall particle at the origin. The other component (\)r\ 
is much more slowly varying and physically incorporates 
the averaged effects of attractive interactions. 

The mean field equation (^) naturally separates into two 
such parts by setting 

<i>m{r) = uo{r) (6) 

and 

(j)Ri{r)=ui{r)+q^,{r), (7) 

with 4>s given by Eq. (^. We call this the basic separation 
and will use it in most of what follows. 

However, other choices can be made. As discussed in 
Appendix B, there exists considerable freedom to vary 
inside the harshly repulsive core region where UQ^r) 
is very large, without affecting the final result when both 
parts of the ERF are taken into account. This flexibility 
can be used to increase the accuracy of approximations 
introduced there that require a slowly varying density re- 
sponse. 

B. First step 

The key idea in the two step method is to compute in an 
initial step the density response to the slowly varying part 
('') of the ERF alone. Physically, this takes account 
of the averaged effects of the attractive interactions mod- 
eled by and we can exploit the fact that the density 
response can be expected to be reasonably slowly varying. 
The response to the remaining harshly repulsive interac- 
tion (jjiiQ is then determined in a second step. 




FIG. 5. Densities for = 0.10 and T = 1.35 for the LJ 
system, the WCA reference system and the nonuniform refer- 
ence system with the ERF determined from Eq. (^) and from 
Eq. (|A1|). The inset focuses in the region around the second 
peak. 

1. Hydrostatic approximation 

Let us first consider the special case where (pj^i varies 
so slowly that it is essentially constant over the range of a 
correlation length in the bulk fluid. The associated density 
Poi^', [4'Ri]ilJ'o) is a functional of the external field (pm and 
a function of the chemical potential and depends only 
on the difference between these quantities. Thus for any 
fixed position ri we can define a shifted chemical potential 

f^o' ^ tJ-o - (I^B.i{ri) (8) 

and shifted field 

^-„\{r)=^ii,{r)-Mr,), (9) 

whose parametric dependence on ri is denoted by a su- 
perscript, and we have for all r the exact relation 

po(r;[0m],Mf)=Po(r;[0Si],AiSO- (10) 

By construction the shifted field (j)]ii{r) vanishes at r = 
ri and it remains very small for r near ri when (j)fji is very 
slowly varying. Thus to determine the density at ri we can 
approximate the R.H.S. of Eq. ( |To| ) by po{ri;[0], ijLq^) = 
Poifv)' density of the uniform fluid (in zero field) 
at the shifted chemical potential /ig^ . We arrive at the 
hydrostatic approximation p8| , p^ for the density arising 
from a very slowly varying field : 

Po{ri;[(j>m],Po) ~ Po(Mo')- (H) 

We refer to po{fiQ^) = p^^ as the tiydrostatic density at ri; 
from Eqs. (||) and it depends only on the /oca/ value of 
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the field (j^m at ri. The hydrostatic approximation is ex- 
act for sufficiently slowly varying c/jjn and has been used in 
more approximate applications of these ideas to hydropho- 
bic interactions in water pof . However in the present ap- 
plication (pfii varies rapidly enough that for quantitative 
accuracy we must use more accurate methods to determine 
the full nonlocal response. 

We now show that the generalized linear response 
method introduced in Ref. ^ provides a simple and ac- 
curate way to determine both the density pm induced in 
the first step as well as the response to 0ijo taken into 
account in the second step. In Appendix B we discuss an 
alternate but somewhat more complicated approach sug- 
gested in Ref. which requires that pm is sufficiently 
slowly varying that gradient type expansions give accurate 
results. 

We start from the exact linear response equation that 
relates small changes in the potential and density for a 
system with external potential 0, chemical potential p, 
and associated density po{r; [(j>],p) = p^{r) [||J|]: 

-W(ri) = j dv2Xo\^u^2M)5p^[r2) (12) 

through the linear response function Xq ^(ri, r2; [p^]) = 
(5(ri-r2)//9^(ri)-co(ri,r2; [p^])- Hereco is the direct cor- 
relation function of the system with density p,f,{r). 



that yields the basic linear response relation (12) for a uni- 
form system. 

This suggests that we can accurately determine the 
desired pm{ri) by using the linear response function 
Xo Po^) of the uniform hydrostatic fiuid in Eq. ( p^ ) 

even when the field (p^j^^ produces significant density 
changes. Assuming a linear density response, we replace 
Sp^ir) in ( p^ ) by the full density change pRi ir) — p^^ , thus 
yielding our final result: 

[pmin) - Po']/Pq' = / dr2Co{ri2;pQ')[pm{r2) - Po'] ■ 



(13) 

Here cq (ri2 ; Pq^ ) is the direct correlation of the uniform 
reference fluid at the hydrostatic density Pq^ . Note that the 
external field appears only implicitly in Eq. (|l^) through 
its local effect on the hydrostatic density Pq^ . 

Equation ( p^ is a linear integral equation relating the 
density pRi{ri) at a given ri on the left side to an integral 
involving the density pm (r2) at all other points and a uni- 
form fluid kernel Ca{ri2 ; Pg^ ) that depends implicitly on ri 
through . This new feature presents no technical diffi- 
culties in determining a numerical solution and Eq. ( p^ ) 
can be solved by any number of standard methods. We 
found that Picard iteration works very well. See Appendix 
C for details. 



2. Linear Response of Hydrostatic Fluid 



C. Second step 



The simple hydrostatic method discussed above approx- 
imates pm{ri) at each ri by the density pg^ of the uniform 
hydrostatic fluid with chemical potential , and thus ig- 
nores the nonlocal effects of the shifted field (j)]^^ on the 
density at ri. To get a more accurate approximation, we 
can use Eq. (^2|) to take into account the linear response of 
the density of the uniform hydrostatic fluid to the shifted 
field . Thus we set p^ = Pq^ and take Scj) = in (|l^). 
This idea was first suggested in Ref. |^ and was shown to 
give accurate results in a number of different applications. 
While a more formal derivation can be given here we 
focus on physical considerations. 

Since (j/j^j (r) is zero at ri , the left side of Eq. ( p^ ) van- 
ishes by construction. We would expect the linear response 
relation between an external field and induced density to 
be most accurate where the field is small — in particu- 
lar where the field vanishes — and at each ri we will use 
the appropriate shifted (hydrostatic) chemical potential 
and shifted field so that this optimal condition continues 
to hold locally. This shift is crucial for the accuracy of 
this method and is its main new feature over previous ap- 
proaches. Moreover, it has been shown that even large 
density fluctuations in a (field free) hard sphere fluid can 
be accurately described using the same Gaussian proba- 
bility distribution that controls small fluctuations j2|] and 



We now determine in a second step the response 
Apfl(r) EE po{r-, [(j)^]) - po{r; [(j)m]) of the relatively slowly 
varying density field po{r: [4>ri]) = PRi{r) to the remain- 
ing harshly repulsive component (pRQ of the ERF, which 
we approximate initially as a hard core of range d. We take 
P4> = P-Ri ill Eq. ( |l^ and again assume a linear density 
response in the "out" region ri > d where the perturb- 
ing potential 0flo(ri) vanishes. This is consistent with the 
simulation results |£1|] showing that the Gaussian prob- 
ability distribution gave a good description even of the 
formation of voids in uniform fluids. 

This linear response assumption gives the approximate 
equation, valid for ri > d: 



0= dr 



f2 Xq 



\ri,r2; [pR,i])l\pR{r2), 



(14) 



where we impose the exact condition po(f2; [07?]) = from 
the hard core interaction for r2 < d in the integration 
over r2. Again we approximate Xo ^ (ri, r2; [p_R,i]) by the 
response function of an appropriately chosen uniform sys- 
tem. As in Sec. 2 and as shown in Ref. H, we find that 



the use of the hydrostatic fluid with density gives ac- 
curate results even when pRi varies rather rapidly. Some 
alternate but less generally useful choices are discussed in 
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Radial distribution functions 50(^2; [fpR 
''^^^])/P^ for two high density states (solid lines with 



FIG. 6. 

= Poir2; [VR 

the origin shifted for clarity) as given by the two step method 
using Eq. (||) for the ERF along with Eqs. (|l|) and com- 
pared to the results of MD simulations (open circles) of the LJ 
fluid. In the inset the two step results are compared directly to 
simulations of the reference fluid in the self-consistent potential 
^MF rpj^g latter tests the accuracy of the two-step method for 
reference fluid correlations induced by the given 4>ji^, while 
the former tests the accuracy of the mean field equation for 
determining 4>^^ ■ 



Appendix B. Thus we arrive at the basic equation for the 
second step of our theory, vahd for ri > d : 



(15) 



Equation ( |l5| ) is a linear equation for Apfl(ri), which 
we can directly solve by iteration or other means. When 
PB.i{f) = , and Co is assumed to vanish for r > d, 
Eq. (15) reduces to the standard PY equation for the 
uniform hard sphere fluid. This has an analytic solution 
IT^J^ and is known to give very good results overall, with 
small errors in the height of the first peak at very high den- 
sities. If still more accuracy is required, we can use modi- 
fied GMSA type equations related to the PY equation 
to describe cq, as discussed in Appendix D. Again we can 
solve Eq. ( p^ ) by iteration. 

This constitutes the second step of our method. The net 
result of this two step process is the desired poir; [4>r]) aris- 
ing from a given 1/)^. This can be substituted into Eq. (|^), 
which can then be iterated to determine the final self- 
consistent and po{r; [0^^^]). See Appendix C for 
further details of the calculations. 




FIG. 7. Radial distribution functions for the lower density 
states indicated compared to the results of MD simulations 
(open circles) of the LJ fluid. The symbols have the same 
meaning as in Fig. 6. 



VI. RESULTS 

We now give a detailed comparison of the radial distri- 
bution functions go(?'2; [0^ ^]) = Po(?'2; ['/'k^D/p^ given 
by the two step method to the results of MD simulations 
of the full LJ fluid. In Fig. (6) we consider two high den- 
sity states with ^ 0.78 and T = 1.35 and p^ = 0.85 
and T = 0.88. At these high densities small errors can be 
seen in the linear response treatment (equivalent to the 
hard core PY equation) of even the uniform hard sphere 
reference fluid. For greater accuracy therefore we used an 
improved GMSA description as briefly described in Ap- 
pendix D. We find by direct comparison with simulations 
of the reference fluid in presence of the self-consistent ERF 
that the two step method indeed gives a very accurate 
description of the nonuniform reference fluid density, as 
illustrated in the inset to Fig. (6). 

The main remaining errors in describing the full LJ fluid 
are thus associated with the mean field approximation for 
the ERF. As already indicated in Fig. (2), this properly 
describes the basic shift in the first peak when compared 
to the WCA uniform reference fiuid, but it also produces 
slight overshoots in the peaks and minima for these high 
density states, as can be seen in Fig. (6). 

Figure (7) gives results for states at T = 1.35 and a 
series of lower densities: p^ = 0.54, p^ — 0.45, and 
p^ = 0.10. Here the linear response treatment of the ref- 
erence system is sufficiently accurate and the attractive 
interactions produce large changes in the correlation func- 
tions. The results seem quite satisfactory, and are compa- 
rable to those given by the best standard integral equation 
methods [0,|lj. Results from the alternate methods dis- 
cussed in Appendix B are equally good and essentially 
indistinguishable on the scale of the graph. 
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VII. FINAL REMARKS 

These results thus give us additional confidence in the 
utility of our general approach. While we certainly do not 
advocate replacing standard and successful integral equa- 
tion methods for the specific problem of the structure of 
uniform simple fluids, these ideas do suggest new ways of 
thinking about some basic issues. We can view the simple 
mean field approximation for the attractive interactions 
along with the generalized linear response treatment of 
correlations in the reference fluid as providing reasonably 
accurate and computationally practical first approxima- 
tions for correlations induced by attractive and repulsive 
interactions. For qualitative and often quantitative work 
they have proved useful in a variety of different applica- 
tions, including cases such as drying near walls where 
attractive forces induce large structural changes and stan- 
dard integral equation methods fail. 

For quantitative accuracy, improvements in both ap- 
proximations may be called for in some cases. Incorpo- 
ration of GMSA type corrections for reference fluid 
correlations is straightforward, as discussed in Appendix 
D, and alternate and probably more accurate treatments 
of the effects of soft cores can be used if needed . Some 
corrections to the simplest mean field equation for the 
ERF, as discussed in Appendix A or in Ref. can also be 
introduced. However, there are some fundamental errors 
arising from the use of any mean field approximation for 
the attractive interactions that cannot be easily avoided. 
The inherent limitations of mean field theory in treating 
long wavelength correlations such as those seen at the criti- 
cal point or arising from capillary waves at the liquid- vapor 
interface are well known. Fortunately in many applica- 
tions of interest such correlations do not play an important 
role, or their effects can be taken into account separately. 
In such cases the ideas discussed here provide a unified 
and physically suggestive perspective capable of giving a 
good qualitative and often a quantitative description of 
the structure of both uniform and nonuniform fluids. 
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APPENDIX A: INTERPOLATED MEAN FIELD 
EQUATION 

Equation is exact as the density p tends to zero, 
where (f>j{^ reduces to the bare field However, the 

next order term [of 0{p)] in a density expansion is incor- 
rect. This can most easily be seen by comparing the known 



density expansions |2) for the LJ system's p{ri; [i^lj]) with 
a LJ particle fixed at the origin, and the reference system's 
Po{ri; [4>r]) with the special wall particle with pair inter- 
action wii(rQi) = 4'R{ri) fixed at the origin. 

We examined an empirical modification of Eq. that 
at low density gives the next term oiO{p) exactly but then 
quickly goes over to Eq. (||) at higher density: 



bT^in) = PMn) - / dr2 {[p,{r,; [c^'i^^]) p^] 



x[l + /(p)/o(ri2)]Fi(ri2;p^)}. 



(Al) 



Here I{p) is an interpolation function that tends to unity 
at low density and to zero at high density, fo{r) = 
exp[— /?uo(^)] ~ 1, and 



Fi(r;p) = [exp{-/3iii(r)/(p)} - l]/I{p). 



(A2) 



Possible choices for I{p) include Ii{p) — {dp/d(5p)o 
= So,p proportional to the reference fiuid isother- 

mal compressibility, and hip) — ^Q p, as suggested by a 
crude argument ||2^ ] based on perturbing the hard sphere 
Ornstein-Zernike equation by a very weak and slowly vary- 
ing potential. At the lowest density studied, p = 0.1, 
Eq. ( |A1| ) with hip) gave a slightly better description of 
the weak second peak than does the simple mean field 
equation (m). See Fig. (5). 



APPENDIX B: ALTERNATE EQUATIONS FOR 
SLOWLY VARYING pm. 

In Sec. we exploited the Gaussian nature of fiuctua- 
tions in the uniform reference system in carrying out both 
steps of the two step method. While this is a good ap- 
proximation for LJ reference system, it may not always 
hold true. In our initial work in Ref. we proposed a 
different, and very general way of carrying out the first 
step, which however requires that pjn varies sufficiently 
slowly that gradient type expansions give good results. 

We started from an exact equation first derived by 
Lovett, Mou and Buff (LMB) ||: 

Vi/9fli(ri) / PBi{ri) = -/3Vi0fli(ri) 

4- J rfr2Co(ri,r2; [pHi])V2Pfli(r2). (Bl) 

The co(ri,r2; [p_Ri]) for a general nonuniform pm is dif- 
ficult to determine, so Eq. (Bl) is generally not very 
useful for practical calculations. However, if pm is rel- 
atively slowly varying, then we can accurately approxi- 
mate co(ri, r2; [pfli]) under the integral in Eq. (Bl) by 
the uniform fluid function co{ri2; P12) p9| ], where pi2 is 
some average density associated with the two points. Then 
Eq. (Bl) can be solved to determine pjn. 

A natural choice for pi2 suggested by a gradient expan- 
sion 111 is pi2 = [pmiri) + pRiir2)]/2 . This gives very 
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good results when pm is reasonably smooth. This is the 
only approximation we make and we can check its accuracy 
by seeing if similar results arise from other approximations 
such as pi2 = PRi{ri) or pi2 = pRi{r2)- Starting with a 
given <j)Ri, we can then solve Eq. (Bl) for the associated 
pRi by iteration, making use of the analytic and accurate 
Percus-Yevick (PY) expressions for the direct correlation 
function of the uniform hard sphere fluid . If more 

accuracy is required we can use GMSA type equations re- 
lated to the PY equation to describe cq . See Appendix 
D. 

In our previous study of the LJ fluid near a hard wall 
0], we used the basic separation of 4'r in Eqs. (||) and 
(13), and found that it indeed produced a slowly varying 
density response. As expected Eq. (Bl) then gave very 
accurate results. However, in the present application, the 
size of the excluded volume region of the fixed particle (of 
order a of the LJ potential) is also the same order as the 
range of the attractive interactions as well as the average 
spacing between particles at high density. If the basic 
separation is used, this "resonance" produces a (j)^^ at 
very high density with small but noticeable oscillations of 
period 27r/cr outside the core and a pronounced minimum 
inside the core at ri = 0. The associated density response 
Poifi'i [0fli]) will have a pronounced maximum at ri = 
and oscillations outside the core, which will cause errors 



in the local expansion method used in Eq. (Bl) and in 
Eq. below. 

Fortunately we can use the flexibility in the choice of 
the field (f>Ri{r) in the core region to produce a much 
smoother density response. More precisely, we can define 
an extended separation of <pR{r) in Eq. (|^) by: 



Mr) -<Po{r) 



and 



c^%Ar)=u^{r)+<i>,{r)+<i>^{r), 



(B2) 



(B3) 



where 4>Q{r) is an essentially arbitrary smooth function 
that is nonzero only in the repulsive core region but with 
(f>o {r) << uo{r), so that (p^Q remains a harshly repulsive 
(essentially hard core) interaction. This separation still 
divides the ERF (j>R{r) into two parts with the physical 
meaning discussed in Sec. but provides some additional 
flexibility in the choice of (f)Ri (r) in the core region that 
can be used to produce a smoother density response in 
the first step. An exact treatment of the response to both 
components of (f>R{r) would of course be independent of 
how the potential was separated. 

Wc found best results by requiring that density response 
to (p^i inside the (hard) core region be constant and con- 
tinuous across the core. In a sense this is the smoothest 
possible choice, at least in the vicinity of the core region. 
This choice can easily be implemented numerically dur- 



iteration and solving for the associated Vi0^^(ri). At 
convergence, the self-consistent PRi{r) is constant inside 
the core and smoother outside the core than that produced 
by the basic separation. 

Using the same extended separation, we have verified by 
comparison with the hydrostatic linear response method 



and with direct simulations that Eq. (Bl) now gives ac- 
curate results for all the states tested here. Thus it offers 
an alternative (though numerically slightly more compli- 
cated) way of carrying out the first step. 

In Ref. [|| we also carried out the second step in a 
slightly different way, effectively combining a local expan- 
sion of pjii with a linear response treatment of the density 
induced by (pRo- In particular, in Eq. ( |l^ we treated 
the (5-function part of Xo ^(ri, i"2; [pfli]) exactly and ap- 
proximated the Co part by the uniform ffuid function at 
the intermediate density pi2 = [pRiiri) + pRi(r2)]/2. The 
accuracy of this approximation can again be checked by 
comparing with other choices such as pi2 = PRiiri). This 
yields the alternate equation for the second step: 

A/9ii(ri)/ofli(ri) = j dv2Co{ri2]pi2)^PR.{r2)- (B4) 

When pRi{ri) is slowly varying, as was the case in all 
the examples studied in the previous work, Eq. (B4) gives 
accurate results, essentially identical to those of Eq. (p^). 
This is also true for most of the states studied in the 
present application, provided that the appropriate ex- 
tended separation is used. In such cases Eqs. (Bl) and 



ing the iterative solution of Eq. (Bl) by simply setting 
Vi/9_Ri(ri) to be zero for all ri inside the core on each 



(B4) can be used as alternate ways of implementing the 
two step method, and for the states shown in Fig. (7) 
they give results on the scale of the graph essentially iden- 
tical to those shown. However for the high density states 
p^ = 0.78 and T = 1.35 and p^ = 0.85 and T = 0.88 in 
Fig. (6) the results using Eq. (B4) vary significantly when 
different choices for pi2 are made. This indicates that for 
these states p^^ is too rapidly varying for Eq. (B4) to 
be trusted. Since Eq. ( p^ ) gives accurate results even for 
these high density states, and makes fewer assumptions 
about the smooth behavior of pRi , it is the preferred way 
to carry out the second step [pO[. 



APPENDIX C: DETAILS OF THE NUMERICAL 
CALCULATIONS 

We give here some details of the numerical solution of 
the basic equations (|), (|l3|), and (H). Eqs. (^) and (g) 
could be used as alternates in the first and second step 
respectively except at the highest densities. We exploit 
the spherical symmetry of the density and the ERF about 
the center of the fixed wall particle, which we take as the 
origin of a spherical coordinate system. Since all these 
equations are used iteratively, we need an efficient and 
accurate method to calculate three dimensional integrals 
over r2 of the general form 
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(CI) 



where due to the spherical symmetry fc(r2) depends only 
on r2 = |r2|, and K{ri2]'p) is a function only of ri2 = |ri — 
Y2 1 and ri and r2 through our choice of the effective density 
p which equals either the hydrostatic density p = or 
the average density p = pi2 = [Pfli(^i) + Pfli(''2)]/2. 
These properties make it advantageous for us to use 



bipolar coordinates [pl[ with the substitution y 



' 12 



2rir2 cos 6' and reduce the three dimensional in- 



tegration to two: 



ri+r2 



I{ri)='^ J dr2r2k{r2) J dyyKiy;p). (C2) 

This transformation is particularly useful if the depen- 
dence of K{y; p) on y is known analytically, since then we 
can explicitly carry out the y integration, and Eq. (02) 
further reduces to a one dimensional integral. All rele- 
vant equations have JsT's that satisfy this condition, thus 
permitting very efficient numerical computations. 
In particular, Eq. (||) becomes: 

n poo 

^R^{r,) - 0i,;(ri) = - / dr2r2 [^0(^2; [4'^]) - P^] 



(C3) 



ri Jo 

ri+r2 



X J dyyui{y), 

In— r2 



while Eq. (na) is transformed to: 



2tt f 

Api?(ri)/poi = — dr2r2ApB{r2) 
ri J 



ri+r2 

X J dyyco{y;Po'). (C4) 

\ri-r2\ 

Equations (|l|) and ( p^ can be similarly rewritten. 

The vector equation (Bl) can be transformed into scalar 
form by taking the scalar product with the unit vector 



r 1 /ri and using the identity 2 ri • r2 
we find: 

dlnpjii{ri) _ 



' 12- 



Thus 



^('^i) 



dn 



dri 

ri+r2 



dpRi{r2) 
dr2 



X / dyy{r^+r2-y-)co{y;p^2)- (C5) 

|i'i-r2| 

In all these equations the integration over the variable 
y can be performed analytically, since we use the PY hard 
sphere direct correlation function Co{y), which is a polyno- 
mial in y 1^, and the ui{y) of the Lennard- Jones potential, 
which is a sum of y~^ and y~^^ terms p^ . Integrals in- 
volving an improved GMSA approximation for co can also 



be carried out analytically. The resulting one dimensional 
integral equations can be solved by Picard iteration, where 
to enforce convergence we use the usual mixing technique. 

In solving these equations the reference potential uq is 
initially taken to be a hard core potential with diameter 
d given by the accurate Verlet- Weiss expressions As 
in the blip function method the result for poir; [(t>R]) 
with Mo approximated by a hard core is linearly extrapo- 
lated into the core region and multiplied by the Boltzmann 
factor of the true soft core potential uq to give the final 
results for po shown in Figs. (6) and (7) . The errors intro- 
duced by this simplified treatment of soft cores are much 
smaller than those arising from our use of the mean field 
approximation for the ERF 0/? . 



APPENDIX D: GENERALIZED MEAN 
SPHERICAL TREATMENT OF REFERENCE 
SYSTEM 

The description of the reference system presented here 
relies on accuracy of the uniform hard sphere fluid direct 
correlation function cq (r) . The generalized linear response 
treatment of Eq. ( p^ ) with p^^ = p^ and cq vanishing out- 
side the core is equivalent to the PY approximation, and 
is surprisingly accurate at intermediate and low densities. 
However, at high density it has noticeable errors, espe- 
cially in the region of the first peak near contact, and for 
quantitative results should be corrected. 

Since the uniform fluid direct correlation function is just 
an input in our approach we can use other, more accurate 
approximations (even results of molecular simulations, if 
such are available). We have found that use of the general- 
ized mean spherical approximation (GMSA) of Waisman 
as implemented by Hoye and StcU |Q, gives consid- 
erable improvement over the original PY approximation. 
Moreover it still preserves the analytic simplicity of the 
resulting co{r) so that the methods of Appendix C can be 
used. 

The GMSA approximates co(r) outside the hard core 
(we set d = 1 here), where PY assumes cq vanishes, by a 
Yukawa function: 



co(r >1) = K- 



(Dl) 



The exact core condition g[r < 1) = then allows one 
to obtain Co(r) inside the core and satisfy the Ornstein- 
Zernike equation: 



co(r < 1) 



^ 3 J 

or r — V- 

2 



, cosh ; 



1 



2A:z2g2 



(D2) 



Requiring consistency between compressibility and virial 
routes for the pressure and agreement with simulations 
gives explicit analytic expressions for a, 6, w, K and z as 
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functions of the packing fraction rj = 7rp/6, as discussed in 
||3^ . We can use this improved expression for co(ri2;Po^) 
in Eq. to describe the hard sphere reference fluid. 

We can also amend our description of the wall particle 
in a similar way by adding the tail correction Eq. (Dl) 
as given by the GMSA to the right side of Eq. (|l|). In 
the absence of attractive forces this equation then reduces 
exactly to a GMSA description of the response of a hard 
sphere fluid to a hard sphere fixed at the origin, and we 
neglect any changes in this correction when attractions are 
taken into account through (j)Ri{r). These GMSA correc- 
tions to the usual linear response treatment are numeri- 
cally significant only for the high density states studied in 
Fig. 6. 
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